arXiv:1509.03321vl [astro-ph.GA] 10 Sep 2015 


Mon. Not. R. Astron. Soc. 000, [JHTG] (2015) Printed 14 September 2015 (MN style file v2.2) 


On the nature of star-forming filaments: II. Sub-filaments 
and velocities 


Rowan J. Smith^*, Simon C. O. Glover^, Ralf. S. Klessen^, Gary A. Fuller^ 

^Jodrell Bank Centre for Astrophysies, Sehool of Physies and Astronomy, University of Manehester, Oxford Road,Manehester, M13 9PL,UK 
Zentrum fiir Astronomie der Univer sit at Heidelberg, Institut fur Theoretisehe Astrophysik, Albert-Ueberle-Str. 2, 69120 Heidelberg, Germany 


14 September 2015 


ABSTRACT 

We show that hydrodynamic turbulent cloud simulations naturally produce large fila¬ 
ments made up of a network of smaller and coherent sub-filaments. Such simulations 
resemble observations of filaments and fibres in nearby molecular clouds. The sub¬ 
filaments are dynamical features formed at the stagnation points of the turbulent ve¬ 
locity field where shocks dissipate the turbulent energy. They are a ubiquitous feature 
of the simulated clouds, which appear from the beginning of the simulation and are 
not formed by gradual fragmentation of larger filaments. Most of the sub-filaments are 
gravitationally sub-critical and do not fragment into cores, however, there is also a sig¬ 
nificant fraction of supercritical sub-filaments which break up into star-forming cores. 
The sub-filaments are coherent along their length, and the residual velocities along 
their spine show that they are subsonically contracting without any ordered rotation 
on scales of ^ 0.1 pc. Accretion flows along the sub-filaments can feed material into 
star forming cores embedded within the network. The overall mass in sub-filaments 
and the number of sub-filaments increases as the cloud evolves. We propose that the 
formation of filaments and sub-filaments is a natural consequence of the turbulent 
cascade in the complex multi-phase interstellar medium. Sub-filaments are formed by 
the high wavenumber, small scale modes in the turbulent velocity field. These are 
then stretched by local shear motions and gathered together by a combination of low 
wavenumber modes and gravitational contraction on larger scales, and by doing so 
build up the extended filaments observed in column density maps. 
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1 INTRODUCTION 


It has become clear that star-forming gas cores are not 
randomly distributed in molecular clouds, but instead are 
threaded along networks of filaments like beads on a string. 
This conclusion be came inescapable fo l lowing the observa¬ 
tions of Herschel (lAndre et aL l2010l: Men’shch ikov et aP 


2010l: Arzouma nian et al.l 12011: iHennemann et al .1 120121: 


■Schneider et al. I I2OI2I : IPalmeirim et al.ll2013lk but had also 
been suggested in previous works (|Schneider_j£^hnng reen 


1979 k and in work with o t her instruments (ISchneider et al 


20K; ljuvela et al.l I2OI2I : iMalinen et al.l 1 20121 : ISeo et al 


2015l k Several authors have shown that filamentary flows 
may play a ke y role in star formation. For instance, 
iKirk et sH (| 20131 ) show that in Serpens South there are mass 
flows of 3.0 X 10“^ M© yr“^ along the filament’s long axis to¬ 
wards the central cluster, and radial mass flows onto the fila¬ 


ment of 1.3 X 10“^ M© yr“^. Similarly. IPeretto et al.l (|2013l ) 
found an inflow rate of 2.5 x 10“^ M© yr“^ along filaments 
towards the central star forming clump in SDC335.579- 
0.272. In numerical simulations, cores embedded within fil¬ 
aments have been shown to for m more massive st ars by ac¬ 
cretion along filamentary flows (|Smith et al.ll201lh . 

Filamentary structure is clearly a key ingredient of 
any theory of star formation, and yet its true nature re- 
mains unexplained . One striking observation is that of 
[Arzoumanian et al.l (| 20 111 ) who propose that filaments have 
a constant width of 0.1 pc in the Herschel Gould Belt sur- 
vey. We add ressed this in the preceding paper in this series, 
ISmith et al.l (2014b), which investigated the projected width 
of filamentary column density features in turbulent clouds. 
However, in this work we will focus on understanding the 
velocity structure of the gas in and around filaments, some¬ 
thing that is crucial for understanding how filaments may 
assemble the mass needed for star formation. 

iHacar et al ] (l2013h showed that although star-forming 
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filaments appear to be continuous features when viewed in 
terms of column density, they actually split up into a net¬ 
work of smaller sub-filaments, or ‘fibres’ when examin ed in 
position-position-velocity space. iTafalla fc Hacarl (|2015l l fur¬ 
ther developed this analysis and showed that these fibres 
were coherent velocity structures with only small velocity os¬ 
cillations along their length. Together they proposed a Tray 
and fragment’ scenario where large filaments were formed 
first, which then fragmented into fibres, with only the dens¬ 
est of these fibres collapsing to form stars. 

In this paper, we show that the formation of a network 
of coherent sub-filaments is a natural consequence of hydro- 
dynamical turbulence. This filamentary network produces 
multiple narrow line compo nents in emission as ob¬ 

served bv iHacar et al.l (|2013h . We define our terminology as 
follows. A filament, or filament network, is a physically re¬ 
lated collection of sub-filaments. A sub-filament is a smaller 
structure, here identified using 3D positions and densities, 
which is part of the larger filament. The connected sub¬ 
filaments will tend to blend together into a single extended 
filament when viewed in column density, as was the case in 
the previous paper in this series. Filamentary structure is 
used to mean any elongated structure that could be either 
a filament or sub-filament. 

The role of filaments in star formation has been 

addressed in previous analyti cal stud ies _ con sidering 

infini t e isothermal cylinders ( Ostriker _ 196jJ: LarsonJ 

Il985l : llnutsuka Mivamal Il992l . 19971 : Iharsonl l2005[ ) and 
i n numerical studies of idealised magnetised cylinders 
(iTillev fc Pudrit3l2003l : iHennebelld l2003|: Seifrie d Wal^ 


(iTillev fc Pudritzll2003l : IHennebelld 120031: ISeifried Walch 
2Q15li and ac cruing cylinders (|Fischera fc Martin 1201? 


Heitschl l2013[ ). However, molecular clouds have a com¬ 


plex density field created by supersonic turbulence 
fe.g. iMa c Low fc Klessen 2004: Elmegree n fc Scald l2004 


IScalo Elmegreenll2004 iKlessen fc Gloveill20l4 ) and rarely 
approach such idealised scenarios. A better approach is 
to simulate clouds as a whole and then consider the fil¬ 
aments arisi ng naturally f ro m su ch turbulent s tructures 
as done in I Jappsen et ^ (l2005ll. ISmith et alJ f20 14bb 


Hennebehe^ Andrd j2013li.[A^eckel fc BurkertI (|20l4 ). and 


Gomez fc Vazouez-Semadenil (|20l4 . We will use this ap¬ 
proach to investigate the filament velocity fields in molecular 
clouds in this paper. 

The paper is structured as follows. In Section [2] we out¬ 
line our numerical techniques and methods. In Section [3] we 
present our analysis of sub-filament velocities, and in Sec¬ 
tion [4] we investigate the time evolution and star formation 
in the filaments. In Section [5] we discuss how our results 
might explain recent observations, and their implication for 
star formation. Finally in Section [6] we give our conclusions. 


2 METHOD 

2.1 Numerical model 


set of mesh-generating points. These mesh-generating points 
and the resulting mesh cells attempt to move with the flow, 
much as the individual particles do in an SPH simulation. 
However, as the mesh is not completely Lagrangian, there is 
generally some residual flux of mass, momentum and energy 
into or out of the cells. These fluxes are computed using a 
Riemann solver, thereby avoiding the need to introduce arti¬ 
ficial viscosity and allowing sharp discontinuities in the flow, 
such as shock fronts, to be modelled with a thickness of only 
1-2 mesh cells. The resulting mesh is adaptable and can be 
refined to give improved resolution in regions of interest. 
This allows the study of problems with an extreme dynamic 
range, that are discontinuous, and that involve fluid insta¬ 
bilities, all while imparting no preferred geometry on the 
problem. 

We model the chemica l evolution of th e gas using the 
NL97 network described in iGlover fc C arki f 2012a), which 
combin es the hydrogen chemistry of Glover fc Mac Lo^ 
(l2007allDl with the simplified treatment of GO form ation 
and destruction introduced in iNelson fc Lang^ (|l997l l. We 
assume that the strength and spectral shape of the ultravi¬ 
olet portion of the interstellar radiation field (ISRF) are the 
same a s the v alues for the solar neighbourhood derived by 
IPrainel (Il978lh note th at this corresponds to a field strength 
of 1.7 iHabind (|l968l l units of 1.6 x 10“^ erg cm“^ s“^. 
We also include the effects of cosmic rays and adopt a rate 
Ch = 3 X 10“^^s“^ for atomic hydrogen, and a rate twice 
the size of this for molecu lar hydrogen. This model was first 
implemented in AREPO in ISmith et al.l (2014a). 

To treat the attenuation of the ISRF due to H 2 self¬ 
shielding, GO self-shielding, the shielding of GO by H 2 , and 
by dust absorption, we us e the treecol algorithm devel¬ 
oped by iGlark et al.l (|2012h . This algorithm computes a Atv 
steradian map of the dust extinction and H 2 and GO col¬ 
umn densities surrounding each AREPO cell, using informa¬ 
tion from the same oct-tree structure that AREPO uses to 
evaluate gravitational interactions between cells. The result¬ 
ing column density map is discretised onto iV pix equal-area 
pixels using the healpix pixelation algorithm (|G6rski et al.l 
I 2 OO 5 I I. In the simulations presented here, we set A^pix = 48. 
To convert from H 2 and GO column densities into the corre- 
spon ding shielding factors, we use shielding functio ns taken 
from lPraine Sz Bertoldil (|l996h and I Lee et al.l (|l996h . respec¬ 
tively. We assume that the radiation field is uniform and 
enters through the sides of the box. The inclusion of time- 
dependent chemistry allows us to calculate the radiative 
heating and cooling of the gas self-consistently within the 
simulation. It is important to model this accurately, as we 
expect filament formation to be significantly easier when 
the effective equation of state of the gas is sub-isothermal 
(|Larsonl [l985l: IPeters et al ]|2 m 3) and so the use of a sim¬ 
ple isothermal or polytropic equation of state may lead to 
unrepresentative results. 


2.2 Simulations 


We per form our simu lations using the moving mesh code 
AREPO E ringe IEmS). This is a quasi-Lagrangian code that 
aims to utilise the strengths of both smoothed particle hy¬ 
drodynamics (SPH) and grid-based adaptive mesh refine¬ 
ment (AMR) codes. The fluid is represented by a series of 
irregular mesh cells that are the Voronoi tessellation of a 


In the previous paper in this series (|Smith et al.l 2014b), 
we considered a set of simulations with different turbu¬ 
lent properties. However, in this paper we chose to fo¬ 
cus on the two simulations containing a natural mix of 
solenoidal and compressive turbulence (e.g. I Schmidt et al.l 
I 2 OO 9 I : iFederrath et ahlboiOl L In both of these simulations. 
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we model 10^ M© gas clouds that have solar metallicity and 
that are initially composed of fully atomic gas. The initial 
condition is that of a uniform sphere with a hydrogen nu¬ 
cleus number density of n ~ 100 cm“^ and a radius of 7 pc. 
This is embedded in a larger 65 pc periodic box contain¬ 
ing a tenuous warm medium with a temperature of several 
thousand Kelvin. The size of this larger box is such that the 
dense cloud at its centre never encounters the box edges. 
This setup is analogous to an isolated molecular cloud in 
the ISM. The turbulence has a P{k) oc k~^ power spectrum 
(specifically, the same spectrum as in runs SI and S2 in the 
nomenclature of Smith et ah, 2014b). The magnitude of the 
root mean square turbulent velocity is normalised such that 
the clouds have an equal amount of kinetic and gravitational 
potential energy at the start of the simulations. Within each 
simulation we focus on the filaments within two different re¬ 
gions, and then follow the time evolution of one of these 
regions over five different snapshots. This gives us twelve 
3.25 pc^ regions to which we can apply our filament finding 
algorithm. 

In order to ensure that the filaments are well resolved 
on all scales we require the Jeans length to be resolved by 
a minimum of 16 cells at all densities. This ensures that 
there is no artificial fragmentation in the gas, and that the 
filament velocities are well resolved at their centres. Where 
necessary, the grid is refined until the condition is satisfied. 
We calculate this Jeans length assuming a fixed temperature 
of 10 K for the gas. Very little of the dense gas is colder than 
this, and most is at least a few K warmer, and so this gives us 
a conservative estimate for the size of the Jeans length. A full 
re solution stud y of the simulations used here was carried out 
in ISmith et al.l (2014b), which also shows the dependence of 
the resolution on the gas density. For example, at densities 
of 10^ cm“^ the cell radius is ~ 5 x 10“^ pc, and at our 
highest densities of 10^ cm“^ the cell radius is 2 x 10“^ pc. 

Star formation is mod elled in the simulation using sink 
particles (iBate et al.l 1995lb T hese were first introduced into 
AREPO by iGreif et al.l (|201lh and we use a slightly modi¬ 
fied version of this routine here. Above number densities of 
10^cm“^, we check whether the densest cell in the simula¬ 
tion and its neighbours satisfy the following three conditions: 
1) the cells are gravitationally bound, 2) they are collapsing 
and 3) the divergence of the accelerations i s less than zero, 
so th e particles will not re-expand (see also iFederrath et al l 
l201Ql b If all these conditions are satisfied, the cell and its 
neighbours are replaced with a sink particle, which interacts 
with the gas cells purely through gravitational forces. Ad¬ 
ditional material can be accreted by the sink particles if it 
is within an accretion radius of race = 0.01 pc (which corre¬ 
sponds to 25 resolution elements at our chosen sink creation 
density) and is bound to the sink. We use relatively large 
sink particles in this study in order to focus on the geom¬ 
etry of the filaments without interference from very dense 
collapsing cores which would distort the average of the fila¬ 
ment profile. Consequently the sinks should not be thought 
of as ‘stars’, but as collapsing cores. 


2.3 Filamentary structure identification 

The filamentary structures we study are identified using the 
DisPer SE (DIScrete P ERsistent Structures Extractor) algo¬ 
rithm (|Sousbiell201lh . This algorithm constructs a Morse- 


Along filament spine 


Vflow 


Vfront 


End on after subtraction 
of COM velocity 



Figure 1. A simple cartoon showing how we decompose the sub¬ 
filament velocities. The flow velocity is the velocity along the spine 
of the sub-filament, and the front velocity is the vector sum of the 
velocity components perpendicular to the spine. The radial and 
rotational velocities are the residual velocities perpendicular to, 
and around the spine after the mean velocity has been subtracted. 


Smale complex from an input density distribution and iden¬ 
tifies the critical points where the density gradient is zero. 
These can then be connected to delineate structures in the 
data. Filamentary structures are found by connecting the 
points such that maxima are connected to saddle-points 
along Morse field lines. To avoid artificial structures being 
identified, we extract only the the structures which have a 
persistence ratio with a probability of 5 sigma or more when 
compared to Poisson noise. We apply DisPerSE to uniform 
3D density grids of regions where filaments were found in 
ISmith et al.l (2014b) and extract the grid points where fila¬ 
mentary structures with a density greater than some thresh¬ 
old occur (2.5 X 10^ cm“^ for the bulk of our cases). We use 
this ‘skeleton’ network to get a series of vectors (hereafter 
called segments) describing the orientation of each filamen¬ 
tary structure. The structures do not have to be linear, but 
instead can bend and twist. For the co-ordinates of the seg¬ 
ments (referred to as spine points), we use the position of 
the densest AREPO cell within a resolution element of the 
DisPerSE input grid for each skeleton point. 

The filament skeleton map obtained above does not con¬ 
tain any information about the properties of the filamentary 
structures, so we have to calculate these from the AREPO 
simulation. For each segment found by DisPerSE, we find 
the mass of the segment by measuring the mass contained 
within a cylinder of radius r = 0.035 pc oriented along the 
vector. This value was chosen for the radius because this was 
found to be the typical flattening radius of a density sub- 
fil ament in our previous paper. To find the flattening radius 
in ISmith et al.l (2014b) we fit a Plummer-like function which 
describes the radial variation in the number density n(r). 
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n(r) = 


Tie 

[1 + 


+ Bsd [cm , 


( 1 ) 


where ric is the central density, i^flat is a radius within 
which the profile is flat, p sets the slope of the power law 
fall-off beyond this radius, and Bsb represents the back¬ 
ground density. This is a common parameterisation of fil¬ 
amentary de nsity profiles and i s widely used in the lit¬ 
erature (e . g iNutter et all l2QQ8l : lArzoumanian et ^ 1201 ll : 
iKirk et ^ l2015l b The typical mean values for the sub¬ 
filaments are ric ~ 10® cm“^, i^flat ~ 0.035 pc, and p ~ 1.9. 
While we c oncentrate on dense star-forming filaments in 
ISmith et al.l (2014b), in the current study we also include 
more diffuse filamentary regions. 

To find the total mass of the filamentary structure we 
add together the mass in each segment within our chosen ra¬ 
dius. Similarly, we add the lengths of the individual segments 
to get the total hlament or sub-filament length. We assign a 
net velocity to each spine point by taking the mass-weighted 
average velocity of all the AREPO cells within this cylinder. 
This means that we are only including the gas at the very 
centre of the sub-filaments where the density profile is flat. 
A different choice of radius, would of course, lead to different 
sub-filament masses and accretio n rates. Howe ver, we chose 
to use the mean value of from ISmith et al.l (2014b) since 
it is a physically motivated choice that is small enough that 
we can be sure only a single sub-filament will be included in 
the mass estimate, and not any nearby neighbour. It should 
be stressed that the total mass of an equivalent filament seen 
in observations will be higher than this value, since it will 
include the total of all the sub-filaments, plus a contribution 
from the more diffuse gas seen in column density out to a 
larger radius. 

In order to investigate the velocity structure of the sub¬ 
filaments, we decompose the velocities along their skeletons 
into the four components shown in Figure [1] The flow veloc¬ 
ity Pflow is the mass-weighted velocity along the sub-filament 
segment, and the front velocity I’front is the mass-weighted 
velocity perpendicular to the sub-filament segment. These 
centre of mass velocity components are then subtracted from 
the velocity in each AREPO cell to find the residual velocities 
around the sub-filament spine. The radial velocity r'rad is the 
mass-weighted velocity measured radially outward from the 
sub-filament centre, as shown in Figure [T] When calculating 
r’rad we include AREPO cells out to a radius of 0.15 pc so 
that gas on the power-law part of the density distribution 
is also included in order to see if mass is flowing towards 
or away from the centre. Finally, the rotational velocity t'rot 
is the mass-weighted velocity component perpendicular to 
Prad, which dcscribcs whether the sub-filament is spinning 
about its long axis. 


3 SUB-FILAMENT VELOCITIES 
3.1 Filament Maps 

Figure [2] shows the column density of four filaments with 
the skeleton obtained from DisPerSE plotted on top.The fil¬ 
aments are named using the convention where S denotes the 
number of the simulation, F denotes the filament within the 
simulation, and T denotes the time in units of 10^ yr since 
the beginning of the simulation. The filament skeletons were 


cut off below a density of 2.5 x 10® cm“® for SIFI, S2FI and 
SIF2, but in the case of S2F2 they were cut off at 2.5 x 10^ 
cm“®. This region was chosen to be larger and more diffuse 
to contrast with the other three cases where some of the sub¬ 
filaments eventually form stars. The sub-filament spines gen¬ 
erally line up well with the column density projection, but 
occasionally a feature is seen in column density that is not 
identified as a sub-filament spine because its density is not 
high enough to meet our selection criteria. For example, a 
long tenuous density feature may have a high column density 
but not have an intrinsic density above our threshold. 

When viewed in another projection the distributions 
look quite different, showing the 3D nature of the filament 
network. If one were to view the column density in the xy 
plane seen in the left panels, then the z direction would be 
the observer’s line of sight. In the panels we see that there 
is the potential for multiple sub-filaments to contribute to 
the emission from a single observed pointing. These overlap- 
pi ng sub-filaments ar e very reminiscent of the observations 
of lHacar et al.l (|20I3l b as we shall discuss further in Section 

El 


3.2 Front and Flow Velocities 

In Figure [3] we show two of the filament skeleton networks, 
but this time colour-code the skeletons so that the most 
massive sub-filaments are red and the least massive are 
dark blue. Arrows show the mean velocity direction and 
relative strength at every tenth spine point. Black crosses 
show where sink particles, representing star-forming cores, 
are formed. Figure [3] clearly shows that the filament net¬ 
works are dynamical features with significant bulk motions. 

A better understanding of how the density sub¬ 
filaments correspond to the turbulent field can be seen in the 
online material q This contains movies showing the velocity 
and density field in slices moving through the simulated fil¬ 
ament networks for SIFI and S2FI. Typically sub-filaments 
are formed at the convergence points of the velocity field 
where mass is swept up by the flow and the density is in¬ 
creased by shocks. This is the formation mechanism of the 
sub-filaments seen in our simulations. In addition to the mo¬ 
tions perpendicular to the sub-filaments (I’front), there are 
also motions parallel to the sub-filaments (i’aow) which can 
be seen particularly clearly in the case of S2FI. 

In Figure |4] we plot the mean front and flow velocities 
for each of the sub-filaments shown in Figure [2l The sub¬ 
filaments are colour coded to show which simulated hlament 
network they are part of. Table [T] shows the mean velocities 
and mass of the sub-hlaments in each hlament network. The 
ratio of the magnitudes of the front velocity to the how ve¬ 
locity is consistently high. By dehnition the front velocity 
must be positive, since it is the magnitude of the two veloc¬ 
ity vectors perpendicular to the hlament long axis, and as 
such we would expect it to be slightly higher, by a factor 
of a/2 , than the how velocity in the case of random sam¬ 
pling. Instead the ratios are much larger than this value, 
about a factor of six, which leads to the conclusion that the 
sub-hlaments are velocity fronts moving through the cloud. 


^ movies are available at http://www.jb.man.ac.uk/~rjs/SI_TI30_g400.mov 
and http://www.jb.man.ac.uk/~rjs/S2_TI70_g400.mov 
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Figure 2. A column density projection of the four different regions within our simulations that we use for our analysis. The column 
density is shown in greyscale and the filament skeleton obtained from DisPerSE is plotted on top. Different sub-filaments are shown in 
different colours so that they can be distinguished. The black boxes show the sub-filaments that are analysed in Figure [6] 
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Figure 3. An illustrative example of the velocities at points along the sub-filament spine for two of our filament networks. The colours 
denote the mass in each sub-filament, with dark blue being the least massive and red the most. The grey arrows show the direction of 
the mass-weighted velocities at points along the sub-filament spines. The black crosses show where sink particles are being formed. 



Front Velocity [km0 

Figure 4. The mean front vs flow velocity for the sub-filaments. 
Black crosses show sub-filaments from SlFl, green asterisks show 
sub-filaments from S2F1, cyan diamonds show sub-filaments from 
S1F2, and pink triangles sub-filaments from S2F2. The sub¬ 
filaments in the grey shaded box are travelling faster than the 
sound speed, which means that there will be shock fronts where 
the sub-filaments collide with slower moving gas. 


The sound speed in the simulated molecular cloud is typi¬ 
cally less that 0.4 km s“^ (about 0.2 kms“^ to 0.3kms“^ in 
the dense gas). Consequently where the sub-filaments move 
through a static or slow-moving gas background they will 
be moving shock fronts. 

If there are velocity gradients along the sub-filament 
spines, there must be a net flow of mass along the sub¬ 
filaments. In Figure [5] we illustrate this using a simple car¬ 



Figure 5. A simple schematic showing how gradients in the ve¬ 
locity along the sub-filaments spine correspond to mass flows. 


toon. Figure [6] shows the variation in the sub-filament ve¬ 
locities along their spine. We plot a sample of connected 
sub-filaments, which are located in the black boxes shown 
in Figure [21 These sub-filaments join smoothly into one an¬ 
other, but are classed as separate objects by DisPerSE as 
there are sharp changes in the angle between their bound¬ 
aries. Each sub-filament is shown in a different colour. The 
lowest panel of the plot shows the density along the filament 
spines. The peaks in the density distribution represent the 
location of dense gas cores. To help us to study the veloci¬ 
ties in the vicinity of these cores, we have highlighted their 
positions with a grey background. 

The top panel in Eigure[6]shows the y-velocity along the 
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Simulation 

^fii 

Mtot 

Mass 

'^front 

1 '^flow 1 

'^rad 

"drot 

crit 

1 Tlfront 1 



[M© ] 

[M© ] 

[kms“l] 

[kms“l] 

[kms“l] 

[kms“l] 


S1F1T140 

47 

194.16 

4.13 (8.83) 

0.75 (0.175) 

0.36 (0.25) 

-0.13 (0.008) 

-0.005 (0.084) 

0.22 (0.23) 

6.22 (11.92) 

S2F1T180 

28 

265.55 

9.48 (9.04) 

0.82 (0.26) 

0.36 5-22) 

-0.11 (0.018) 

0.014 (0.071) 

1.04 5-70) 

5.39 53-74) 

S1F2T140 

38 

196.17 

5.16 (7.05) 

1.25 5-39) 

0.49 5-42) 

-0.11 (0.008) 

0.012 5-088) 

0.37 5-36) 

7.12 5i-52) 

S2F2T100 

35 

63.18 

1.81 5-18) 

1.29 (0.67) 

0.37 (0.42) 

-0.07 5-011) 

0.003 5-069) 

0.03 5-02) 

6.75 5i-50) 


Table 1. The mean value and its standard deviation (in parenthesis) of the mass and velocities of all the sub-filaments in the simulated 
regions. is the number of sub-filaments, Mtot the total mass in sub-filaments in the network, and crit is the critical line mass ratio 
in Equation [ 3 ] The mean absolute magnitude of Uflow is shown so that it can be fairly compared to Ufront- The ratio | Ufront/'^flow I is 
found by first calculating the ratio for each individual sub-filament and then taking the average. 



Figure 6. The velocity along one line-of-sight, Uflow? '^^rad and density along one axis for a series of connected sub-filaments. Each 
sub-filament is shown by a different line colour. The grey shading highlight the positions of peaks in the density profile along the spine, 
which correspond to the locations of dense cores, and the thin grey line shows the trend of the velocities. There are oscillations in the 
velocity field along the filament spine, but generally there is both radial contraction and a flow of material along the sub-filaments onto 
the cores. 


sub-filaments from S1F1T303, and the x-velocity for sub¬ 
filaments in S2F1T390. We chose to plot x- or ^-velocities 
depending on the orientation of the box in Figure O The 
second panel shows the decomposed component Ufiow along 
the filament spine. In both cases there are oscillations along 
the spines (the two cases are not identical because the sub¬ 
filaments are not completely straight but have kinks and 
twists). S uch oscillatory motion s have been previously ob¬ 
served bv iHacar fc Tafall^ (|201lh where they correlated well 
with the locations of dense cores. In subsequent studies, how¬ 


ever, oscillations were still seen, but they n o longer corre¬ 
lated with the cores (jTafalla Hacarll2015h . In our simu¬ 
lated sub-filaments the oscillations and cores are not well 
correlated, implying that they are transient motions that 
vary depending on the random turbulent field in which the 
sub-filaments are embedded. We also calculated the Fourier 
transform of Ufiow along the filament spines and found no 
dominant frequency mode. Again, this suggests that the 
oscillations are driven by the surrounding turbulent field. 
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Sub-filament Mass [M,„|] 


Figure 7. The net mass flow along individual sub-fllament spines. 
Black crosses show sub-filaments in SlFl, green asterisks show 
sub-filaments in S2F1, cyan diamonds show sub-filaments in 
S1F2, and pink triangles show sub-filaments in S2F2. 


rather than being a consequence of a periodic gravitational 
instability along the filament spine. 

In addition to the oscillations, there are also larger scale 
trends in the velocity field, shown by the thin grey lines 
in Figure [6l For example, in the sub-filament in S1F1T303 
shown in red, there is a net positive gradient along the sub¬ 
filament. By examining the cartoon shown in Figure [5] we 
can see that this corresponds to a net flow of mass along the 
sub-filament towards a core at the end. The adjoining blue 
sub-filament has a negative net gradient along its length, 
which corresponds to a mass flow towards the centre. It is 
therefore clear that mass can be transferred along the sub¬ 
filaments. The cores tend to appear preferentially where the 
large-scale velocity field changes and where there is a flow 
of mass towards the core. This demonstrates the importance 
of gravity in forming the cores. Depending on the relative 
strength of the tidal field compared to the self-gravity of the 
sub-filament, the sub-filament can be interpreted either as 
an accretion flow that feeds an existing over-density at the 
end of the sub-filament, or instead it may fragment to form 
its own cores. We will discuss sub-filament fragmentation 
further in Section [H 

To estimate the magnitude of such flows we cannot just 
sum the mass moving along the sub-filaments, because if we 
were to take the scenario in Figure [5] the net total would 
come to zero in both cases despite there being strong mass 
flows. Instead we fit a straight line to the flow velocity along 
every sub-filament spine weighted by the mass in each seg¬ 
ment for each of the four simulations shown in Figure [2] The 
gradient of this line is then used to obtain a crude order-of- 
magnitude estimate of the net flow of mass per parsec using 
the relationship Maow = where Mai is the total 

mass of the filament. Figure [7| shows the estimated value of 
Maow for each sub-filament in the four simulations, plotted 
against its mass. 



Figure 10. The mean radial vs rotational velocity for all the 
sub-filaments. Black crosses show sub-filaments from SlFl, green 
asterisks show sub-filaments from S2F1, cyan diamonds show sub- 
filaments from S1F2, and pink triangles show sub-filaments from 
S2F2. 


The mass flows are both negative and positive, mean¬ 
ing that in some sub-filaments mass is flowing towards 
the centre, and in others to the end. For the dense sub¬ 
filaments, tens of solar masses of gas can be moved along the 
sub-filament per Myr. For comparis on in the super-cr itical, 
southern filament of Serpens South, iKirk et ^ (|2013h mea¬ 
sured a gas flow of 30 M© Myr“^ along the filament, which 
is of the same order of magnitude as the values we mea¬ 
sure in our simulations. Large-scale filaments have also been 
observed to have conver ging motions attributed to collapse 
along their long axis fe.g IZernickel et al.ll2013l : IPeretto et al.l 
12014 ). However, it is important to stress that this process 
can be both constructive and destructive, as sub-filaments 
with a positive Maow will become increasingly diffuse at 
their centre. The arrow of time does not point inevitably 
towards increasing densities, and many sub-filaments will 
never form stars. Instead they may either be sheared apart 
by turbulence, or cannibalised by nearby accreting cores. 

3.3 Radial and Rotational Velocities 

Having investigated the net velocities along the sub-filament 
spines, we now turn our attention to the local residual ve¬ 
locities. The radial and rotational velocities were obtained 
by subtracting off the centre-of-mass velocity for each sub¬ 
filament segment and then decomposing the velocities as dis¬ 
cussed in Section [2.31 In Figure [6] we also plot the mean ra¬ 
dial velocity for each segment along the sub-filament spines. 
In almost all cases the radial velocity is between 0 and -0.2 
kms“^ indicating subsonic radial contraction of gas onto the 
spine. The filaments are, therefore, accreting radially as well 
as along the filament spine. As before, the radial contraction 
is not uniform, but instead varies along the sub-filament. 

To understand the details of the residual motions we 
show in Figures [8] and [9] the radial and rotational velocities 
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Figure 8. Illustrative slices showing the radial velocity in some sub-filament cross-sections. Red points show the position of cells with a 
net outward velocity and blue shows cells moving inwards. 
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Figure 9. Illustrative slices showing the rotational velocity in some sub-filament cross-sections. Red points show the position of cells 
with an anti-clockwise velocity and blue shows cells moving clockwise. 


from three sub-filament segments in SlFl and S2F2 as il¬ 
lustrative examples. The projection is face-on so that the 
spine vector points into the page. Figure [8] shows the radial 
velocity of the gas cells in the segments. In each case there 
is a high velocity flow of the more diffuse material onto the 
sub-filament spine centre, but at th e centre the velocities 
are more quiescent. As discussed in ISmith et al.l (2014b), 
the sub-filament sections are elongated, not circular, which 
shows that the sub-filaments are not true cylinders. The flow 
of gas onto the sub-filaments is more reminiscent of a col¬ 
liding flow than of a radial inflow and there is still some gas 
moving away from the sub-filament centre. Figure [9] shows 
the rotational velocity of the gas cells in the segments. There 
is no strong evidence for uniform rotation about the spine 
centre. While there are velocities perpendicular to the radial 
direction these are not ordered at scales of ~ 0.1 pc (roughly 
the sonic scale) and seem to be from random turbulence or 
shearing in the field rather than rotation. 

Figure fTOl shows the mass-weighted mean radial and ro¬ 
tational velocities of all the sub-filaments. For each simu¬ 
lation the radial velocities tend to cluster around the same 


value, since the sub-filaments are embedded in the same tur¬ 
bulent bulk flow. For both Urad and Urot, the velocities are 
subsonic despite being measured over a r = 0.15 pc region. 
This is due to most of the mass being concentrated at the 
sub-filament centre where it can be supported by thermal 
motions (a volume-weighted mean does indeed give a large r 
velocity). A similar result was found in ISmith et al.l (|2012l T 
where the velocity dispersion measured in simulated N 2 H^ 
emission from cores embedded in filaments was narrower 
than that calculated by a volume-weighted average. We will 
discuss the observational consequences of the coherency of 
the sub-filaments velocities in Section[5]and in future papers. 
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Figure 11. The evolution of SlFl and S2F2 at intervals of 2.17 x 10^ years (ten snapshot intervals) as seen in the x-y plane. The colours denote the mass in each sub-filament, with 
dark blue being the least massive and red the most massive. Sink particles are shown by crosses. The grey arrows show the direction of the filament velocity at every tenth filament 
spine point. While the overall shape of the network stays similar the number and positions of the sub-filaments do change over time. 
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4 TIME EVOLUTION AND 

FRAGMENTATION 

The filamentary networks presented in the last section were 
studied at a single point in time. However, in order to inves¬ 
tigate how the sub-filaments are involved in star formation, 
we need to consider their dynamical evolution. Figure ITT] 
shows the evolution of the filament network in SlFl and 
S2F1. While the overall outline of the network stays similar, 
the number of sub-filaments and their positions vary. 

In Table [2] we show how the masses of the sub-filaments 
and the number of sink particles evolve. The time quoted for 
each snapshot is measured from the beginning of the simula¬ 
tion. The number of sub-filaments shows no consistent trend 
with time, which emphasises that filamentary substructure 
is a general feature of molecular clouds. The mass in the 
sub-filaments, on the other hand, does clearly increase as 
the cloud evolves. This is because the mass in the region 
is increasing due to global converging motions. The con¬ 
verging motions were initially the result of the large-scale 
turbulent modes gathering gas together, but over time these 
regions become gravitationally unstable, causing the gas to 
collapse and the density to increase. Over the studied pe¬ 
riod the mass in the sub-filament spines in SlFl increased 
by 2.15xlO“^Mo yr^^ and in S2F1 by 9.7xlO“^Mo Myr^^ 
Substantial amounts of mass can therefore be gathered into 
dense coherent sub-filaments over the lifetime of the cloud, 
making it more susceptible to fragmentation. 

The mass per unit length is a useful way of quantifying 
whether a filament can remain in a stable equilibrium. Us¬ 
ing a model of an infinite self-gravitating cyl inder with no 
magnetic support. Ilnutsuka Sz Mivamal (|l997l i found that a 
filament cannot be supported thermally below the critical 
line mass 

Merit = 2c^/G = 16.7 Mo pc-\ (2) 

where Cs is the sound speed, G is the gravitational constant 
and T is the temperature. We calculate the critical ratio for 
our filaments using the relation 

Crit — Tfline/Tfcrit 5 (3) 

where Mune is the measured mass per pc along the filament 
spine. 

As a consequence of the sub-filament mass increasing, 
the mean critical ratio of the sub-filaments rises as well, 
and so does the number of super-critical sub-filaments. In 
all of the simulated regions, most of the sub-filaments re¬ 
main sub-critical and do not fragment. However, multiple 
sink particles (representing star-forming cores) are formed 
in the sub-filaments that are super-critical. Th i s in a gree- 
ment with the observations of iTafalla Hac^ (|2015l i who 
find that star formation only occurs in the dense ‘fertile fi¬ 
bres’ seen in N 2 H^ and not in any of the sub-filaments only 
detected in C^^O. We will discuss this further in Section [5Tl 

One could get a different estimate of the critical line 
mass if the internal velocity dispersion were to be added to 
the sound speed. However, our definition of the sub-filament 
mass is a very conservative one, as we only include material 
in the flat part of the density profile in Equation [T] to avoid 
double counting mass where the sub-filaments bend or are in 
close proximity. Within this region the residual velocities are 
subsonic as discussed in Section [T3l and so do not contribute 


Simulation 

'^front 

1 '^flow 1 

^rad 

'drot 


[kms“^] 

[ km s“^] 

[kms“^] 

[kms“^] 

S1F1T238 

0.74 

0.38 

-0.100 

-0.030 

S1F1T260 

0.70 

0.33 

-0.089 

0.008 

S1F1T281 

0.71 

0.30 

-0.102 

0.010 

S1F1T303 

0.75 

0.36 

-0.128 

-0.005 

S1F1T325 

0.89 

0.46 

-0.119 

0.014 

S2F1T346 

0.82 

0.33 

-0.122 

-0.012 

S2F1T368 

0.78 

0.31 

-0.132 

-0.019 

S2F1T390 

0.82 

0.36 

-0.118 

0.014 

S2F1T411 

0.68 

0.41 

-0.086 

0.017 

S2F1T433 

0.63 

0.35 

-0.083 

0.005 


Table 3. The evolution of the mean velocity components as 
the simulations evolve. The mean absolute magnitude of Uflow is 
shown so that it can be fairly compared to Ufront • 


significant support. Furthermore, as shown by Figures [6] and 
M the radial velocity component is preferentially directed 
inwards making it unclear how much support such motions 
offer against collapse. In the simulations we see that cores 
(i.e. sink particles) are only formed in super-critical sub¬ 
filaments, which suggests that our definition is reasonable. 

In Section [3] we discussed the decomposed velocities in 
the sub-filaments. Table [3] shows that there is no systematic 
change in the velocities as the sub-filaments evolve, with 
the exception of rotation. As shown in Figure [9l there is no 
clear evidence for ordered rotation in the filament, and so 
the variation in Urot probably comes from contributions of 
random turbulent motions. The other velocity components 
appear to be set by global properties of the cloud, explaining 
the lack of variation. The front and flow velocities are de¬ 
termined by the large-scale velocity gradients within which 
the filament network is embedded, and the low value of Urad 
is because on scales corresponding to the width of the sub- 
filaments, we expect the velocity dispersion to be subsonic 
(|Hever fc Bruntlt2004l b 


5 DISCUSSION 

5.1 Comparison to observations 

So far we have concentrated on a purely theoretical analysis 
of the sub-filament velocities, but it is also important to de¬ 
termine whether our results are consistent with observations. 
We will do this more fully in future work where we plan to 
do a detailed comparison of the derived filament velocities 
in the observational plane. However, this paper would be 
incomplete if we did not show that our sub-filaments actu¬ 
ally resemble real position-position-velocity structures seen 
in molecular clouds. 

Our appro a ch is largely motivated by the observations 
of iHacar et al ] ll2013h who show that the Taurus filament 
exhibited multiple narrow line components when observed in 
emission. It is this finding that allows them to deduce 
the existence of velocity coherent sub-filaments or ‘fibres’ 
within the clouds. 

We have shown that hydrodynamic turbulent clouds 
contain a network of sub-filaments in agreement with ob¬ 
servations, but we have not yet proven that this leads to 
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Filament 

Time 

A^fii 

Mean mass 

Total mass 

Mean crit 

Kcrit 

Max crit 

A^sink 


[Myr] 


[M© ] 

[M© ] 





S1F1T238 

2.38 

30 

3.82 (3.90) 

114.65 

0.20 (0.14) 

0 

0.62 

0 

S1F1T260 

2.60 

31 

4.48 (8.40) 

138.86 

0.18 (0.14) 

0 

69 

0 

S1F1T281 

2.81 

61 

3.06 (6.70) 

186.52 

0.22 (0.21) 

0 

0.90 

0 

S1F1T303 

3.03 

47 

4.13 (8.83) 

194.16 

0.22 (0.23) 

3 

1.46 

1 

S1F1T325 

3.25 

56 

5.39 (8.47) 

301.87 

0.57 (1.34) 

7 

9.48 

2 

S2F1T160 

3.46 

33 

5.77 (8.25) 

190.47 

0.38 (0.52) 

2 

2.58 

1 

S2F1T170 

3.68 

34 

6.58 ©28) 

223.71 

0.53 (0.62) 

5 

3.65 

6 

S2F1T180 

3.90 

28 

9.48 (9.05) 

265.56 

1.04 (1.70) 

8 

7.18 

9 

S2F1T190 

4.11 

20 

13.01 (15.27) 

260.24 

1.20 (1.72) 

6 

5.54 

16 

S2F1T200 

4.33 

26 

10.61 (15.32) 

275.80 

1.54 (3.99) 

7 

20.1 

17 


Table 2. The evolution of the sub-filament properties for SlFl and S2F1. Nf[i is the number of sub-filaments, crit is the critical line 
mass in Equation O A^crit is the number of sub-filaments with a critical ratio greater than one, Max crit is the greatest value of crit for 
the sub-filaments, and A^sink is the number of sink particles in the filament. The standard deviation in each mean is shown in brackets 
to give a feel for the variation in the sample. 


the observed emission profiles. To test this hypothesis, we 
perform radiative transfer calculations of S1F1T281 and 
S2F1T368 using the radmc-3d cod43- RADMC-3 d models 
a variety of radiative processes including dust thermal emis¬ 
sion and absorption, dust scattering and atomic and molec¬ 
ular line emission. It is this latter capability which is of key 
use here. We use our AREPO simulations as input for the den¬ 
sity, temperature, velocities, and CO abundance, and then 
use RADMC-3D to calculate the level populations and per¬ 
form a ray-trace. 

As is not necessarily in local thermodynamic equi¬ 

librium (LTE) throughout our filamen ts, we use the l arge ve¬ 
locity gradient (LVG) approximation (|Sobolevlll957l h which 
assumes that gas which is widely separated in position is also 
widely separated in velocity space, meaning that only local 
radiation will affect the level populations jj As RADMC-3D 
cannot work directly with the unstructured Voronoi mesh 
used within AREPO, we first map the simulation data onto 
a 400^ Cartesian grid with a pixel size of 0.008 pc before 
producing our synthetic emission maps. The C^^O isotope 
is not tracked directly in the chemical model, and so we as¬ 
sume that it h as an abundan ce that is 557 times smaller than 
that of ^^CO (|Wilsonlll999l i. To make the final image and 
spectra we then smooth the image with a gaussian beam. To 
test the effect of beam size we use three different beams with 
full width half maxima of 2, 4, and 8 pixels respectively. At 
the distance of Taurus (140 pc), these would correspond to 
beams of size ~ 0.4', 0.8' and 1.6'. 

We do not adopt a preferred distance for our model, but 
instead calculate the intensity per pixel in physical units, not 
the angular size. This can be converted to a pixel scale in 
degrees once a choice is made as to the distance of the object. 
Since we are not modelling a specific region here we use the 
more general case, so that the reader can scale the values 
to the distance of their preferred object. In order to make 
the lines more easily comparable to brightness temperature, 
we convert our intensities into Kelvin by multiplying by the 

^ This code is publicly available at the website www.ita.uni- 
heidelberg.de / ~dullemond / software / radmc-3d / 

^ For a full description of how the LVG app r oxima tion has been 
implemented in RADMC-3D, see lShettv et al 


factor A^/2fc, where A is the wavelength of the line and k is 
Boltzmann’s constant. 

Figure [12] shows the resulting integrated emission map 
for the J = 1 ^ 0 transition of C^^O, based on simulation 
S1F1T281 viewed in the x-y plane (i.e. the observer’s sight¬ 
line is the z-axis in Fig[2|) with a 4 pixel FWHM beam. For 
illustrative purposes we show the line profiles at two points 
along the filament seen in the integrated emission. The line 
profiles show multiple velocity components in the C^^O (1-0) 
line, in agreement with observations. Moreover, the line com¬ 
ponents are narrow due to the individual sub-filaments being 
coherent with subsonic velocities in the dense gas along the 
sub-filament spine as discussed in Section [3.31 

Figure [13] explores the origin of these lines. The bottom 
panel shows the density along the observer’s line of sight 
and the top shows the corresponding velocity, for the two 
line profiles in Figure [12] The lines are colour-coded with 
their velocity to make it easy to see the velocity of the high 
density peaks along the sight line. We find there are at least 
as many density peaks as there are visible velocity com¬ 
ponents, meaning that the number of observed components 
gives a lower limit on the amount of substructure in a region 
similar to the one that we simulate. 

We plan to carry out a full analysis of the peak veloci¬ 
ties and line-widths that would be observed in our simulated 
clouds in future work where we will treat the problem in de¬ 
tail and include observational effects such as noise. However, 
in order to get an idea of how common multiple line com¬ 
ponents might be, we perform the following analysis. Using 
the synthetic position-position-velocity (PPV) cube gener¬ 
ated by RADMC-3D, we can obtain a line profile for each 
position in the image. When the line profile is bright (here 
defined as the peak emission being above 2 K), we find the 
number of maxima in the profile. This brightness threshold 
is arbitrarily chosen to ensure that we are only including pix¬ 
els which lie in the filament identified in the column density 
map. 

Table [4] shows the percentage of bright pixels with mul¬ 
tiple maxima in their line profiles for different simulations, 
line transitions and beam sizes. Some general trends emerge. 
The percentage of lines with multiple maxima varies between 
the simulated clouds, which is to be expected as they have 
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Figure 12. Synthetic emission map for the (1-0) line, for simulation SlFl observed with a 4 pixel beam. The panels on 

the right show the line profiles from the two highlighted regions. In agreement with observations there are multiple narrow velocity 
components seen in the line. 
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Figure 13. The density and line-of-sight velocity along the beam for the two line profiles shown in Figure [TJ] The colour scale also 
shows the line-of-sight velocity to make it easier to see the velocity at which the density field peaks. Typically the number of components 
in the line gives a lower limit to the number of over-densities in the beam. 


different filament networks. There are more lines with mul¬ 
tiple components in the lower energy J = 1 ^ 0 line than 
the higher energy J — 2 ^ 1 line. This is because some 
of the sub-critical sub-filaments have densities much lower 
than the J — 2 ^ 1 critical density, and so the lines from 
these sub-filaments are faint. The number of line compo¬ 
nents also decreases as the beam size increases because the 
lines become broader and are more likely to be blended to¬ 


gether. Regardless of this, the presence of multiple maxima 
in the emission lines is a ubiquitous feature of the synthetic 
observations. 

It is not only in the predicted line emi ssion that the fil¬ 
ament s are in agreement with observations. iTafalla HacaJ 
Jiolil) showed that core formation was only occurring in a 
subset of ‘fertile fibres’ and that many of their observed fi¬ 
bres were not forming stars. When a fibre did form a core. 
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pation of turbulent modes with short wavelengths and that 
the filaments seen in column density are formed by long 
wavelength turbulent modes gathering the gas together. 


Filament 

Transition 

Beam 

[pix] 

2 maxima 

[%] 

More than 2 

[%] 

S1F1T281 

1-0 

2 

34.12 

11.18 

S1F1T281 

1-0 

4 

32.82 

7.35 

S1F1T281 

1-0 

8 

28.77 

3.80 

S1F1T281 

2-1 

2 

31.90 

8.13 

S1F1T281 

2-1 

4 

28.16 

4.70 

S1F1T281 

2-1 

8 

21.41 

1.18 

S2F1T368 

1-0 

2 

25.29 

2.59 

S2F1T368 

1-0 

4 

19.63 

1.43 

S2F1T368 

1-0 

8 

12.72 

0.38 

S2F1T368 

2-1 

2 

21.06 

2.14 

S2F1T368 

2-1 

4 

16.36 

0.95 

S2F1T368 

2-1 

8 

9.80 

0.41 


Table 4. The percentage of line profiles in the simulated 

PPV cube with a peak brightness greater than 2 K that have 
multiple maxima. 


it tended to form as part of a chain. In the simulations of 
collapsing turbulent clouds we observe similar behaviour. 
Most of the sub-filaments are sub-critical, but when they do 
reach criticality multiple sink particles (representing collaps- 
in g cores) tend to form, in a greement with the predictions 
of llnutsuka Mivamal (|l997l b Only very few sub-filaments 
become super-critical within 10® years of their formation. 
Most others need considerably longer to accrete sufficient 
mass - if they manage at all. This fact may help to account 
for the low efficiency of star formation in molecular clouds. 

5.2 Fray and gather 

From their observations iTafalla Hacarl (|2015l l proposed a 
Tray and fragment’ scenario where a filamentary density en¬ 
hancement was first formed by two gas flows colliding, which 
then Trayed’ into sub-filaments, and then fragmented into 
cores. In our simulations we have the advantage of being 
able to follow the time evolution of the gas. In the collapsing 
turbulent clouds the Tray’ stage of the process where sub¬ 
filaments are formed occurs first This can be seen from the 
fact that the number of sub-hlaments can either increase or 
decrease as the filament evolves (see Table [2|), and that even 
in diffuse gas such as S2F2 there are many sub-hlaments. 
Dense gas has this morphological appearance right from the 
beginning. These sub-hlaments are then gathered together 
by large-scale motions within the cloud. Initially these are 
the large-scale turbulent modes, which have bulk velocities 
on the scale of the box length. However, as the turbulence 
decays the cloud becomes more and more gravitationally un¬ 
stable which also gather s the sub-hla ments together. This 
process was illustrated in ISmith et al.l (2014b) when we dis¬ 
cussed the evolution of S2F1. 

The formation of sub-hlaments ca n also be thought 
of in terms of the t urbulent cascade ([Kolmogorovl Il94ll : 
iBoldvrev et al.l l2002l l . I Schmidt et al.l (|2008h analysed high- 
resolution turbulent boxes with various forcing scales and 
found that the most dissipative structures are intermediate 
between hlamentary and sheet-like structures. We propose 
that the sub-hlaments are formed from the turbulent dissi¬ 


5.3 Caveats and future work 

This work has so far neglected the effects of magnetic helds. 
This was a deliberate choice in order to understand the 
velocity he l d ari sing from turbulence and gravity alone. 
iKirk et ^ (|2015l l studied the density prohles of hlaments 
formed in simulations of MHD turbulence and found that 
magnetic helds provided additional pressure support which 
broadened the column density prohles and re duced frag¬ 
menta tion compared to the hydrodynamic case. iHennebelld 
(|2013h studied MHD simulations of non-self-gravitating hl¬ 
aments and suggested that ion-neutral friction may play a 
key role in determining the width of hlaments. However at 
present there is no investigation of possible differences be¬ 
tween hlament velocities and sub-structure in the hydrody- 
namical and magnetohydrodynamical cases. This is some¬ 
thing that we hope to address in future work. 

We have also mainly conhned ourselves to a theoret¬ 
ical analysis, with the exception of Section 15.11 where we 
tested whether the simulations were compatible with obser¬ 
vations. In future work we plan to carry out a fuller anal¬ 
ysis, where we will test how synthetic observations of these 
clouds correspond to the simulations and how we may best 
derive the physical properties of the clouds from the ob¬ 
servations. In previous work on hlaments, we showed that 
star-forming cores embedded within hlaments would have 
very different observational properties f r om th ose expected 
for isolated spherical cores. ISmith et al.l (|2012h showed that 
the turbulent velocities in the diffuse envelopes of the hla¬ 
ments could obscure the blue asymmetr ic line prohle signa¬ 
ture expected from a collapsing core. In ISmith et al.l (|2013h 
we demonstrated that for massive protostars forming at the 
centre of a collapsing hlamentary hub, the optically thick 
line prohles lacked self-absorption features and there were 
multiple optically thin line components in dense gas tracers. 
In simple spherical collapse models one would expect strong 
self-absorption and a single narrow optically thin compo¬ 
nent at the line centre. We expect that mock observations of 
the sub-hlament networks simulated in this paper will once 
again reveal line emission that could not have been predicted 
from simple analytical models of isolated star formation. 


6 CONCLUSIONS 

In this paper we have used high resolution AREPO simula¬ 
tions to show that turbulence in collapsing clouds naturally 
leads to the formation of a network of sub-hlaments in agree¬ 
ment with recent observations of Taurus. Our main conclu¬ 
sions are as follows. 

(i) Filaments seen in column density maps actually con¬ 
sist of a network of sub-hlaments when viewed in 3D. 

(ii) The hlaments are dynamical features that move ac¬ 
cording to the large-scale turbulent velocity held in the 
molecular cloud. These bulk velocities are predominantly 
perpendicular to the long axis of the hlament. 
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(iii) There are flows of mass along the spine of the sub¬ 
filaments of the order of tens of solar masses per Myr. These 
can compress the sub-filament, increasing its mass and mak¬ 
ing it more likely to fragment. Alternatively the sub-filament 
may be sheared apart by turbulence, or form an accretion 
flow which feeds mass into a nearby star-forming region. 

(iv) At the dense spine of the filament, where the den¬ 
sity profile is flat, the gas velocities are both subsonic and 
coherent. The radial velocities are of order —0.1 kms“^ di¬ 
rected inwards, and there is no ordered rotation about the 
sub-filament spines. 

(v) Most s ub-filaments are sub - critica l according to the 
formalism of llnutsuka fc Mivamal (|l997[ ) and do not frag¬ 
ment into cores. In the sub-filaments which are super¬ 
critical, multiple cores are forming, in agreement with ob¬ 
servations. 

(vi) Sub-filaments are a general feature of turbulent 
molecular clouds. They appear at all times and for all gas 
densities. 

(vii) The mass in sub-filaments and the number of super¬ 
critical sub-filaments increases as the molecular clouds 
evolve. In our simulated clouds the mass in sub-filament 
spines increased at a rate of 215 M© Myr“^ in SIFI, and by 
97 M© Myr“^ in S2FI. 

(viii) As the network of sub-filaments overlap along a typ¬ 
ical observed sightline, this results in multiple narrow line 
components appearing in simulate d line emissio n, in 

agreement with the observations of lHacar et al.l (|20I3l i. 

(ix) We propose that the formation of filaments and sub¬ 
filaments is a natural consequence of a turbulent cascade. 
Sub-filaments are formed by the high wavenumber, short 
wavelength modes. These are then gathered together by 
large-scale, low wavenumber modes and gravitational col¬ 
lapse to form the extended filaments observed in column 
density maps. 
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